#include <string>
#include <deque>
#include <algorithm>
#include <iostream>
#include <vector>
#include <sstream>
#include <boost/lexical_cast.hpp>
#include <algorithm>
#include "dynamic_matrix.h"
#include "stlhelper.h"
#include "portable_timer.h"
#include "progress_indicator.h"
#include "pg_copy_guard.h"
#include "pg_db_cursor.h"
#include "index_string.h"
//# define DEBUG_LEVEL 19
//# define DEBUG_LEVEL 19
#include "debug_func.h"
#include "get_arg.h"
using std::cerr;
using std::cerr;
using std::deque;
using std::string;
using std::istream;
using std::cout;
using std::cin;
using std::swap;


//________________________________________________________________________________________

// t_sparse_matrix

//		only used matrix positions are filled and flagged with matrix_filled

//________________________________________________________________________________________
struct t_sparse_matrix
{
	t_dynamic_matrix		matrix;
	t_dynamic_matrix_bool	matrix_filled;
	void set(unsigned i, unsigned j, double val)
	{
		matrix.set(i, j, val);
		matrix_filled.set(i, j, true);
	}
	bool less_than(unsigned i, unsigned j, double val) const
	{
		return matrix_filled.at(i, j) && matrix.at(i, j) < val;
	}
	t_sparse_matrix(): matrix(0.0), matrix_filled(false){}
	bool get_element(unsigned i, unsigned j, double& val) const
	{
		if (matrix_filled.at(i, j))
		{
			val =  matrix.at(i, j);
			return true;
		}
		return false;
	}
};




void print_cluster(const deque<unsigned>& clusters,
					const t_sparse_matrix& sparse_matrix, const t_index_string& gene_names)
{
	for (unsigned i = 0; i < clusters.size() - 1; ++i)
	{
		for (unsigned j = i + 1 ; j < clusters.size(); ++j)
		{
			double val;
			if (sparse_matrix.get_element(i, j, val))
			{
				cout	<< gene_names.string_at(i)		<< "\t"
						<< gene_names.string_at(j)		<< "\t"
						<< val							<< "\n";
			}
		}
	}
}


//________________________________________________________________________________________

// gather extra cluster members

//________________________________________________________________________________________
void gather_extra_cluster_members(	t_program_args&					args,
									deque<deque<unsigned> >& 		clusters,
									const t_sparse_matrix&			sparse_matrix,
									const t_index_string&			gene_names,
									double							ks_cutoff,
									unsigned&						cnt_clustered_pairs,
									unsigned&						cnt_clustered_genes)
{
	deque<deque<unsigned> > 		cluster_new_members(clusters.size());
	unsigned sz = gene_names.size();
	t_progress_indicator	dots (args.VERBOSE(), 50);
	//iterate through clusters
	for (unsigned clust = 0; clust < clusters.size(); ++clust)
	{
		// iterate through cluster members
		for (unsigned member = 0; member < clusters[clust].size(); ++member)
		{
			// clustered gene id
			unsigned gene = clusters[clust][member];

			dots++;

			// now test all in matrix which have relationships with this clustered gene
			for (unsigned test = 0; test < sz; ++test)
			{
				// no self comparisons
				if (test == gene)
					continue;
				// stored in i < j order (triangle not square matrix) to save space
				unsigned ii = gene;
				unsigned jj = test;
				if (ii > jj)
					swap(ii, jj);
				if (sparse_matrix.less_than(ii, jj, ks_cutoff))
					cluster_new_members[clust].push_back(test);
			}
		}
	}
	
	deque<unsigned> clustered_genes;
	// add new clusters to old and uniquify
	for (unsigned clust = 0; clust < clusters.size(); ++clust)
	{
		// save all genes which have been clustered at all
		clustered_genes.insert(	clustered_genes.end(),
								cluster_new_members[clust].begin(),
								cluster_new_members[clust].end());
		clusters[clust].insert(	clusters[clust].end(),
								cluster_new_members[clust].begin(),
								cluster_new_members[clust].end());

		sort(clusters[clust].begin(),  clusters[clust].end());
		remove_duplicates(clusters[clust]);
		for (unsigned i = 0; i < clusters[clust] .size(); ++i)
			cout	<< clust << "\t"
					<< gene_names.string_at(clusters[clust][i]) << "\n";
	}
	
	// get the number of clustered genes
	sort(clustered_genes.begin(),  clustered_genes.end());
	remove_duplicates(clustered_genes);
	cnt_clustered_genes = clustered_genes.size();

	for (unsigned clust = 0; clust < clusters.size(); ++clust)
	{
		cout << ">" << clust << "\n";
		print_cluster(clusters[clust], sparse_matrix, gene_names);
	}


}

//________________________________________________________________________________________

// retrieve_ks_pairs

//		only used matrix positions are filled and flagged with matrix_filled

//________________________________________________________________________________________
bool retrieve_ks_pairs	(	t_pg_cursor&		pg_db,
						 	t_program_args&		args,
							const string& 		ks_table,
							t_sparse_matrix&	sparse_matrix,
							t_index_string&		gene_names_index,
							unsigned& 			cnt_pairs)
{
	debug_func_cerr( 8 );
	string sql_cmd = "SELECT "
							"gene_id1, gene_id2, ks "
					 "FROM "
							+ ks_table;
	

	
	
	// create binary cursor
	if (!pg_db.Declare(sql_cmd, true))
		return pg_print_sql_err("create a cursor to retrieve value and gene pairs",
								pg_db,
								sql_cmd, args.CERR());
	
	{
		t_progress_indicator	dots (args.VERBOSE(), 30000);
		unsigned batch_size = 5000;
		// get data in batches
		while (pg_db.Fetch(batch_size))
		{
			int curr_batch_size = pg_db.Tuples();
			if (!curr_batch_size)
				break;
			for (int i = 0; i < curr_batch_size; ++i)
			{
				++cnt_pairs;
				++dots;
				
				double val = pg_db.getvalue<double>(i, 2);
				
				// save ks pair
				string gene_id1 = pg_db.getvalue(i, 0);
				string gene_id2 = pg_db.getvalue(i, 1);
				unsigned gi1 = gene_names_index.index(gene_id1);
				unsigned gi2 = gene_names_index.index(gene_id2);
				if (gi1 > gi2)
					std::swap(gi1,  gi2);
				sparse_matrix.set(gi1, gi2, val);
			 }
	
		}
	
	}
	
	// close cursor
	if (!pg_db.Close())
		return pg_print_sql_err("close the  cursor for "
								"create a cursor to retrieve ks and gene pairs",
								pg_db, "CLOSE CURSOR", args.CERR());
	
	return true;
	
}



//________________________________________________________________________________________
void retrieve_clusters(	deque<deque<unsigned> >& clusters, t_index_string& gene_names, istream& is)
{
	while (1)
	{
		unsigned cluster;
		string member;
		is >> cluster >> member;
		if (!is)
			break;
		bool added = false;
		unsigned member_id = gene_names.index(member, added);
		if (added)
		{
			cerr	<< "cluster " << cluster <<
						" contains an unrecognised gene name ["
					<<	member << "]\n";
		}
		
		while (clusters.size() <= cluster)
			clusters.push_back(deque<unsigned>());
		clusters[cluster].push_back(member_id);
	}
	
	for (unsigned i = 0; i < clusters.size(); ++i)
		if (!clusters[i].size())
			cerr << "Error\n\tCluster " << i << " has no members!\n";
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_portable_timer timer;
	
	// don't use C streams
	std::ios::sync_with_stdio(false);

	// saving parameters
	string exe_parameters;
	for (int i = 0; i < argc; ++i)
		exe_parameters += string(argv[i]) + ' ';

	// we use cout to print progress so flush per call
	// cout.setf(std::ios_base::unitbuf);
    char buff[65536];
    cout.rdbuf()->pubsetbuf(buff, 65536); // enable buffering

	// get command line arguments
	t_program_args args;
	if (!GetArguments(argc, argv, args))
		return 1;

	
	//____________________________________________________________________________________
	// 		Connect to PostgreSQL
	//____________________________________________________________________________________
	t_pg_cursor pg_db(args.db_login.c_str(), "panda_cursor");	//"dbname=all_vs_all_gene_bl2seq user=postgres"
	if ( pg_db.ConnectionBad() )
	{
		cerr	<< "\nERROR!! Connection to database ["
							<< args.db_login << "] failed.\n"
							<< pg_db.ErrorMessage() << "\n\n";
		return 1;
	}
	
	
	//____________________________________________________________________________________
	// 		get all pairs
	//____________________________________________________________________________________
	t_index_string gene_names(500000);
	t_sparse_matrix			sparse_matrix;
	unsigned cnt_pairs = 0;
	args.VERBOSE() << "\tRetrieving Ks pairs...\n";
	if (!retrieve_ks_pairs		(	pg_db, args, args.ks_table,
									sparse_matrix, gene_names, cnt_pairs))
		return 1;
	unsigned cnt_genes				= gene_names.size();
	args.VERBOSE() << "\t" << cnt_pairs << " pairs retrieved.\n";
	args.VERBOSE() << "\t"<< cnt_genes << " genes retrieved.\n";

	//____________________________________________________________________________________
	// 		retrieve clusters
	//____________________________________________________________________________________
	args.VERBOSE() << "\tRetrieving current cluster members...\n";
	deque<deque<unsigned> > clusters;
	retrieve_clusters		(clusters, gene_names, cin);
	args.VERBOSE() << "\t" << clusters.size() << " clusters retrieved.\n";
	
	unsigned cnt_clustered_pairs	= 0;
	unsigned cnt_clustered_genes	= 0;
	args.VERBOSE() << "\tGather extra cluster members...\n";
	gather_extra_cluster_members(args, clusters, sparse_matrix, gene_names, args.ks_cutoff,
								cnt_clustered_pairs, cnt_clustered_genes);

	
	args.VERBOSE() << string(80, '_') << "\n";
	args.VERBOSE() << cnt_clustered_pairs << " pairs out of " << cnt_pairs << " could be clustered\n";
	args.VERBOSE() << cnt_clustered_genes << " genes out of " << cnt_genes << " could be clustered\n";
	args.VERBOSE() << timer.pretty_reset() << " Seconds\n";

	return 0;
}






